function Tree=df_setup(N,P,MAT)
  
  L = round(log(N/P)/log(2))+1;
  h = 1/N;
  
  %1. process the sparse matrix
  [ii,jj] = ndgrid(0:N-1);
  T = zeros(N,N); %center
  edge = find(mod(ii,P)==0 | mod(jj,P)==0);
  T(edge) = 1;
  vert = find(mod(ii,P)==0 & mod(jj,P)==0);
  T(vert) = 2;
  
  [I,J,S] = find(MAT);
  tmp = find(T(I)>=1 & T(J)>=1 & T(I)+T(J)<4);  S(tmp) = S(tmp)/2;
  tmp = find(T(I)>=2 & T(J)>=2);  S(tmp) = S(tmp)/4;
  NEW = sparse(I,J,S); %matrix weighted corrected
  
  GRD = reshape(1:N*N,N,N);
  Tree = cell(L,1);
  %---
  if(1)
    ell = 1;
    W = 2^(ell-1)*P;
    nck = N/W;
    Tree{ell} = cell(nck,nck);
    for i=1:nck
      for j=1:nck
        if(i<nck)
          is = (i-1)*W+[1:W+1]; ia = is(1); ib = is(end);
        else
          is = [1 (i-1)*W+[1:W]]; ia = is(2); ib = is(1);
        end
        if(j<nck)
          js = (j-1)*W+[1:W+1]; ja = js(1); jb = js(end);
        else
          js = [1 (j-1)*W+[1:W]]; ja = js(2); jb = js(1);
        end
        %maps
        glbidx = GRD(is,js);
        [ii,jj] = ndgrid(is,js);
        gud = ones(size(ii)); %all good
        lclidx = zeros(W+1,W+1);
        lclidx(find(gud==1)) = 1:numel(find(gud==1));
        isin = (ii~=ia & ii~=ib & jj~=ja & jj~=jb);
        LINS = lclidx(find(gud==1 & isin==1));
        GINS = glbidx(find(gud==1 & isin==1));
        LBDS = lclidx(find(gud==1 & isin==0));
        GBDS = glbidx(find(gud==1 & isin==0));
        %construct mat
        %ALL = NEW(glbidx,glbidx);
        %invA = inv(full(ALL(LINS,LINS)));
        %B = ALL(LBDS,LINS);
        %S = full(ALL(LBDS,LBDS)) - full(B*invA*B');        if(issparse(S)==1) error('wrong'); end;
        ALL = NEW(glbidx,glbidx);
        invA = inv(full(ALL(LINS,LINS)));        invA = (invA+invA')/2;
        B = ALL(LBDS,LINS);
        S = ALL(LBDS,LBDS) - B*invA*B';        S = (S+S')/2;
        Tree{ell}{i,j} = struct('LINS', LINS, 'LBDS', LBDS, 'GINS', GINS, 'GBDS', GBDS, 'invA', invA, 'B', B, 'S', S); clear invA B S;
      end
    end
  end
  %---
  for ell=2:L-1
    W = 2^(ell-1) * P;    Q = W/2;
    nck = N/W;
    Tree{ell} = cell(nck,nck);
    for i=1:nck
      for j=1:nck
        if(i<nck)
          is = (i-1)*W+[1:W+1]; ia = is(1); ib = is(end); im = ia+Q;
        else
          is = [1 (i-1)*W+[1:W]]; ia = is(2); ib = is(1); im = ia+Q;
        end
        if(j<nck)
          js = (j-1)*W+[1:W+1]; ja = js(1); jb = js(end); jm = ja+Q;
        else
          js = [1 (j-1)*W+[1:W]]; ja = js(2); jb = js(1); jm = ja+Q;
        end
        %maps
        glbidx = GRD(is,js);
        [ii,jj] = ndgrid(is,js);
        gud = (ii==ia | ii==im | ii==ib | jj==ja | jj==jm | jj==jb);
        lclidx = zeros(W+1,W+1);
        lclidx(find(gud==1)) = 1:numel(find(gud==1));
        isin = (ii~=ia & ii~=ib & jj~=ja & jj~=jb);
        LINS = lclidx(find(gud==1 & isin==1));
        GINS = glbidx(find(gud==1 & isin==1));
        LBDS = lclidx(find(gud==1 & isin==0));
        GBDS = glbidx(find(gud==1 & isin==0));
        %construct mat
        nb = numel(find(gud==1));
        ALL = zeros(nb,nb);
        cut = ( (ia<=ii&ii<=im) & (ja<=jj&jj<=jm) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i-1,2*j-1}.S; Tree{ell-1}{2*i-1,2*j-1}.S=[];
        cut = ( (im<=ii|ii==ib) & (ja<=jj&jj<=jm) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i  ,2*j-1}.S; Tree{ell-1}{2*i  ,2*j-1}.S=[];
        cut = ( (ia<=ii&ii<=im) & (jm<=jj|jj==jb) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i-1,2*j  }.S; Tree{ell-1}{2*i-1,2*j  }.S=[];
        cut = ( (im<=ii|ii==ib) & (jm<=jj|jj==jb) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i  ,2*j  }.S; Tree{ell-1}{2*i  ,2*j  }.S=[];
        %Schur complement
        invA = inv(ALL(LINS,LINS));        invA = (invA+invA')/2;
        B = ALL(LBDS,LINS);
        S = ALL(LBDS,LBDS) - B*invA*B.';        S = (S+S')/2;
        Tree{ell}{i,j} = struct('LINS', LINS, 'LBDS', LBDS, 'GINS', GINS, 'GBDS', GBDS, 'invA', invA, 'B', B, 'S', S); clear invA B S;
      end
    end
  end
  %---
  if(1)
    ell = L;
    W = 2^(ell-1)*P;    Q = W/2;
    nck = N/W;    if(nck~=1) error('wrong'); end
    if(1)
      if(1)
        i=1;
        j=1;
        is = 1:W;    ia = 1;    ib = 1;    im = 1+Q;
        js = 1:W;    ja = 1;    ib = 1;    jm = 1+Q;
        glbidx = GRD(is,js);
        [ii,jj] = ndgrid(is,js);
        gud = (ii==ia | ii==im | jj==ja | jj==jm);
        lclidx = zeros(W,W);
        lclidx(find(gud==1)) = 1:numel(find(gud==1));
        isin = (ii>ia & jj>ja);
        LINS = lclidx(find(gud==1 & isin==1));
        GINS = glbidx(find(gud==1 & isin==1));
        LBDS = lclidx(find(gud==1 & isin==0));
        GBDS = glbidx(find(gud==1 & isin==0));
        %construct mat
        nb = numel(find(gud==1));
        ALL = zeros(nb,nb);
        cut = ( (ia<=ii&ii<=im) & (ja<=jj&jj<=jm) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i-1,2*j-1}.S; Tree{ell-1}{2*i-1,2*j-1}.S=[];
        cut = ( (im<=ii|ii==ib) & (ja<=jj&jj<=jm) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i  ,2*j-1}.S; Tree{ell-1}{2*i  ,2*j-1}.S=[];
        cut = ( (ia<=ii&ii<=im) & (jm<=jj|jj==jb) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i-1,2*j  }.S; Tree{ell-1}{2*i-1,2*j  }.S=[];
        cut = ( (im<=ii|ii==ib) & (jm<=jj|jj==jb) );        idx = lclidx(find(gud==1 & cut==1));
        ALL(idx,idx) = ALL(idx,idx) + Tree{ell-1}{2*i  ,2*j  }.S; Tree{ell-1}{2*i  ,2*j  }.S=[];
        %Schur complment
        invA = inv(ALL(LINS,LINS));        invA = (invA+invA')/2;
        B = ALL(LBDS,LINS);
        S = ALL(LBDS,LBDS) - B*invA*B.';        S = (S+S')/2;
        S = inv(S);        S = (S+S')/2;
        Tree{ell}{i,j} = struct('LINS', LINS, 'LBDS', LBDS, 'GINS', GINS, 'GBDS', GBDS, 'invA', invA, 'B', B, 'S', S); clear invA B S;
      end
    end
  end
  